## Base RF model
set.seed(42)
library(ggplot2) # for autoplot() generic
library(dplyr)
library(sf)
library(stars)
library(tmap)
dat_sf <- st_read("./data/glacier_clim.shp")
## Reading layer `glacier_clim' from data source `/Users/u0784726/Dropbox/Data/devtools/glacier_mb/data/glacier_clim.shp' using driver `ESRI Shapefile'
## Simple feature collection with 95086 features and 31 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 67.47845 ymin: 27.4918 xmax: 103.8915 ymax: 45.35089
## Geodetic CRS: GCS_unknown
Load and crop basemap
basemap <- read_stars("~/Dropbox/Data/summer/natural_earth/HYP_50M_SR_W/HYP_50M_SR_W.tif")
basemap2 <- st_crop(basemap, st_bbox(dat_sf))
basemap2 <- st_rgb(basemap2)
projstr <- "+proj=lcc +lon_0=83.0 +lat_0=32.0 +lat_1=32.0 +lat_2=38.0 +datum=WGS84 +no_defs"
t2t2m_2008 <- read_stars("~/Dropbox/Data/summer/harr/output/t2_2008_amjjas.nc")
st_crs(t2m_2008) <- projstr
t2m_2018 <- read_stars("~/Dropbox/Data/summer/harr/output/t2_2018_amjjas.nc")
st_crs(t2m_2018) <- projstr
Make delta
t2m_d <- t2m_2018 - t2m_2008
t2m_trend <- read_stars("~/Dropbox/Data/summer/harr/output/t2_b_ann.nc")
st_crs(t2m_trend) <- projstr
tptp_2008 <- read_stars("~/Dropbox/Data/summer/harr/output/prcp_2008_ann.nc")
st_crs(tp_2008) <- projstr
tp_2018 <- read_stars("~/Dropbox/Data/summer/harr/output/prcp_2018_ann.nc")
st_crs(tp_2018) <- projstr
Make delta
tp_d <- tp_2018 - tp_2008
Seasonality
tp_seas <- read_stars("~/Dropbox/Data/summer/harr/output/prcp_2018_seas.nc")
st_crs(tp_seas) <- projstr
t2mx <- c(t2m_2008, t2m_2018)
names(x) <- "2m Air Temp (AMJJAS)"
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(x) +
tm_raster(palette = "-inferno",
alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
Plot \(\delta\)
names(t2m_d) <- "2m Air Temp (AMJJAS) (change 2018 - 2008)"
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(t2m_d) +
tm_raster(alpha = 0.75,
style = "cont") +
tm_layout(legend.outside = TRUE)
print(m1)
names(t2m_d) <- "2m Air Temp (AMJJAS) (change 2018 - 2008)"
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(t2m_d) +
tm_raster(alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
Plot \(\delta\) with mass balance (continuous scale)
m1 <- tm_shape(t2m_d) +
tm_raster(alpha = 0.75,
style = "cont") +
tm_shape(dat_sf) +
tm_symbols(col = "mb_mwea",
size = 0.25,
alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
names(t2m_trend) <- "2m Air Temp Trend (slope 1991-2020)"
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(t2m_trend) +
tm_raster(palette = "-PuOr",
alpha = 0.75,
style = "cont", midpoint = 0) +
tm_layout(legend.outside = TRUE)
print(m1)
## stars_proxy object shown at 1093 by 537 cells.
tpx <- c(tp_2008, tp_2018)
names(x) <- "Total PPT (Annual)"
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(x) +
tm_raster(palette = "YlGnBu",
alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
Plot \(\delta\)
names(tp_d) <- "Total PPT (Annual) (change 2018 - 2008)"
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(tp_d) +
tm_raster(alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
Plot \(\delta\) with mass balance
m1 <- tm_shape(tp_d) +
tm_raster(alpha = 0.75,
style = "quantile") +
tm_shape(dat_sf) +
tm_symbols(col = "mb_mwea",
size = 0.25,
alpha = 0.75,
style = "quantile") +
tm_layout(legend.outside = TRUE)
print(m1)
tp seasonalitynames(tp_seas) <- "PPT seasonality (JJA / ANN)"
m1 <- tm_shape(basemap2) +
tm_raster(palette = "Greys", alpha = 0.5) +
tm_shape(tp_seas) +
tm_raster(palette = "-PuOr",
alpha = 0.75,
style = "quantile", midpoint = 0) +
tm_layout(legend.outside = TRUE)
print(m1)
## stars_proxy object shown at 1093 by 537 cells.